Fermions on bilayer graphene: symmetry breaking for B = and v = 0. 
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We extend previous analyses of fermions on a honeycomb bilayer lattice via weak-coupling renor- 
malization group (RG) methods with extremely short-range and extremely long-range interactions 
to the case of finite-range interactions. In particular, we consider different types of interactions 
including screened Coulomb interactions, much like those produced by a point charge placed either 
above a single infinite conducting plate or exactly halfway between two parallel infinite conducting 
plates. Our considerations are motivated by the fact that, in some recent experiments on bilayer 
grapheneii^ there is a single gate while in others^, there are two gates, which can function as the 
conducting planes and which, we argue, can lead to distinct broken symmetry phases. We map out 
the phases that the system enters as a function of the range of the interaction. We discover that the 
system enters an antiferromagnetic phase for short ranges of the interaction and a nematic phase 
at long ranges, in agreement with Refs. @ and @. While the antiferromagnetic phase results in a 
gap in the spectrum, the nematic phase is gapless, splitting the quadratic degeneracy points into 
two Dirac cones each^. We also consider the effects of an applied magnetic field on the system 
in the antiferromagnetic phase via variational mean field theory. At low fields, we find that the 
antiferromagnetic order parameter, A(B) — A(0) ~ B 2 . At higher fields, when lo c > 2Ao, we find 
that A(B) ~ w c /[ln(a; c /A(0)) + C], where C ~ 0.67 and uj c = eB/m*c. We also determine the 
energy gap for creating electron-hole excitations in the system, and, at high fields, we find it to be 
auj c + 2A(B), where a is a non-universal, interaction- dependent, constant. 



I. INTRODUCTION 

The problem of interacting fermions on an A-B stacked 
honeycomb bilayer lattice has been of great interest, both 
theoretically and experimentally, due in no small part to 
its band structure. Typically, the band structure pos- 
sesses two inequivalcnt points, labeled K and K' = K, 
in the Brillouin zone at which the two low-energy bands 
make contact at four Dirac-like cones^. However, in 
the case where we neglect all but nearest-neighbor hop- 
ping terms, these four Dirac cones merge into a single 
quadratic degeneracy point£. The density of states at the 
Fermi level at half-filling therefore becomes finite in this 
special case. This leads to logarithmic divergences in the 
uniform and static susceptibilities to various symmetry- 
breaking orders at zero temperature^, and thus we ex- 
pect such phases to appear for arbitrarily weak interac- 
tions over some finite range of temperatures. 

One key motivation for the study of the honeycomb 
bilayer in particular is the fact that bilayer graphene, 
a material of great interest, possesses this lattice struc- 
ture. In particular, measurements of the conductivity of 
suspended bilayer graphene in various external electric 
and magnetic fields as well as different carrier densities^, 
followed by compressibility measurements^, provide ev- 
idence suggesting that a symmetry-breaking state ap- 
pears; the two possible states argued for in these works 
are a quantum anomalous Hall phase, in which the sys- 
tem develops a non-zero quantized Hall conductivity in 
the absence of a magnetic field, and a nematic phase, in 
which the quadratic degeneracy points each split into two 
massless Dirac-like cones. A more recent experiment^ 
finds evidence for the presence of a nematic phase^ 



by measuring the width of a peak in the resistivity as a 
function of the carrier density at different temperatures 
and by measuring cyclotron gaps as a function of the 
applied magnetic field for different filling factors. An- 
other, even more recent, experiment^ uses measurements 
of two-terminal conductance to argue for a state in which 
the system develops an energy gap in its spectrum, in ap- 
parent disagreement with the previous experiment, since 
the nematic state is gapless. Two other experiments^^ 
also conclude that a gap due to symmetry breaking is 
present, but are inconclusive about the exact nature of 
the state. Yet another experiment^ performs measure- 
ments on multiple samples, finding a bimodal distribu- 
tion of both conducting and insulating samples. 

In addition to experimental investigations, there have 
also been a number of theoretical studies regarding 
the nature of the symmetry-breaking states of the A- 
B stacked honeycomb bilayer. One such investigation-^ 
uses variational methods to argue for a ferromagnetic 
phase for long-range interactions and a calculation of 
the susceptibility towards an antiferromagnetic phase 
within an RPA approximation, followed by a mean-field 
calculation of the associated order parameter, to ar- 
gue for said phase for short-range interactions. Two 
other works employ a mean-field approach to argue for 
a "(layer) pseudospin magnet" stated and a ferromag- 
netic stated, while later investigations argue for a layer- 
polarized stated within the mean-field approximation 
and a quantum anomalous Hall stated with mean field 
plus Gaussian fluctuations. Another work by the same 
group^ considers the phase of the system within the 
mean-field approximation as a function of external elec- 
tric and magnetic fields and finds a quantum anomalous 



Hall state, a quantum Hall ferromagnetic insulator, and a 
layer-polarized state. Another worki^ uses Hartree-Fock 
methods to argue for a layer-polarized state in the ab- 
sence of a magnetic field, and argues for the existence of 
an anomalous quantum Hall state in the presence of a 
magnetic field. A very recent investigation by the same 
group arrives at a "(layer) pseudospin antiferromagnetic" 
state for the honeycomb bilayer—. Another very recent 
investigation finds, using Hartree-Fock methods, a coex- 
istence of a quantum spin Hall state and a layer-polarized 
state that may be turned into a pure layer-polarized 
state with the application of a sufficiently strong electric 
field^i. The work of Ref. [22| treats the problem of the 
combined effect of an in-plane magnetic field, i.e., with a 
Zeeman term only, and an applied perpendicular electric 
field, within a self-consistent mean field approximation. 

The mean-field approach, as a means of predicting the 
low-temperature phase of the system, however, has one 
major disadvantage. It does not treat the leading loga- 
rithmic divergences that appear in perturbation theory 
correctly, and therefore may not lead to results that are 
correct, even qualitatively^. On the other hand, the 
renormalization group (RG) approach does. An example 
of this point is provided by the treatment of interacting 
fermions on a one-dimensional chain^. In the spinless 
case (see, for example, Ref. HH), mean field theory pre- 
dicts that there will be a charge density wave state for 
arbitrarily weak interactions at half filling, while RG pre- 
dicts a Luttinger liquid for weak interactions. We there- 
fore believe that RG is more accurate as a means of pre- 
dicting the state that our system enters than a mean field 
approach. There are now several papers that employ the 
RG approach in studying the honeycomb bilayer in the 
weak-coupling limit£~— 

A paper by one of us and Yang£ uses this method to 
argue for the existence of a nematic state for extremely 
long-range interactions in the case of spin-^ fermions, 
and a similar paper by Lemonik et. aZ.— ^ which fol- 
lowed comes to the same conclusion. Another paper by 
one of us£ investigates the extremely short-range limit, 
namely the repulsive Hubbard model, using the same ap- 
proach, and finds that a system of spin-i fermions will 
enter an antiferromagnetic state; it is also argued within 
that this state should persist even in the strong-coupling 
limit. This has been very recently confirmcdSS using a 
combination of quantum Monte Carlo and functional RG 
methods. 

In the present work, we will extend the analyses 
conducted in Refs. i and i to the case of finite- 
ranged density-density interactions in the case of spin-^ 
fermions, and to finite temperature, following the meth- 
ods described in detail in Ref. H3. We are interested in 
determining how the system transitions from the anti- 
ferromagnetic state to the nematic state as we increase 
the range of the interaction. In particular, we consider 
two types of interactions, namely 1) screened Coulomb 
interactions, much like those produced by a point charge 
placed exactly halfway between two parallel infinite con- 



ducting plates, 2) same as 1) but for a single infinite 
conducting plate. Each form considered includes a pa- 
rameter £ that is proportional to the range of the inter- 
action. 

Our procedure for the calculation is as follows. We 
start with a tight-binding lattice model for the honey- 
comb bilayer with a density-density interaction between 
electrons. From this, we may derive a low-energy effec- 
tive field theory, which takes the same form as that given 
in the previous work&; said field theory contains 9 dif- 
ferent contact quartic couplings that arc allowed by the 
symmetries of the honeycomb bilayer lattice and by Fierz 
identities. We find the values of the coupling constants 
in terms of the interaction in our original tight-binding 
lattice microscopic model. We then use these values as 
the initial conditions for the RG flow equations derived 
to one-loop order in Ref. l27l which we integrate numer- 
ically. We then perform the analysis outlined therein to 
determine the phase that the system enters. 

For very short-ranged interactions, we find an antifer- 
romagnetic phase, in agreement with the previous work&. 
This is determined by monitoring the susceptibilities to- 
ward various phases as the temperature is lowered; in this 
case, the only susceptibility that diverges at the critical 
temperature is towards the antiferromagnetic phase. If 
we increase the range, we will enter a region where the 
susceptibilities toward both the antiferromagnetic and 
nematic phases both diverge, although not necessarily 
with the same exponential. This transition occurs for 
values of £ anywhere from about 0.4 to 2 lattice spac- 
ings, depending on the form and overall strength of the 
interaction. Finally, upon increasing £ further, we enter 
a pure nematic phase, again in agreement with the previ- 
ous work££. This happens for values of £ anywhere from 
4 to about 10 lattice spacings, again depending on the 
form and overall strength of the interaction. 

For the dipole-likc potential, much like the one that is 
produced by a point charge a distance £ above a single 
infinite conducting plate, 
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where Uq is a sufficiently small constant, so that our sys- 
tem is at weak coupling. For the potential produced by a 
point charge placed exactly half-way between two infinite 
conducting plates separated by a distance 2£, 



V(v) = 4U A "c 
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where K n {x) is a modified Bcssel function of the second 
kind and the last equation holds asymptotically at large 
values of r/£. Note that these expressions, as they are, 
diverge at r = — that is, they diverge for two electrons 
on the same lattice site. Because of this, and due to the 
fact that the electrons in the tight-binding model cannot, 
strictly speaking, be considered point particles, we ren- 



der these on-site interactions finite by setting said inter- 
actions equal to a constant A times the nearest-neighbor 
interaction; our main results were found for A = 1.2; this 
choice makes our interactions monotonically decreasing 
with distance. 

In each case, we find a qualitatively identical phase 
diagram, except that the values of £ at which the transi- 
tions described earlier occur are smaller for the one-plate 
case than for the two-plate case. 

Another way to think of our results is as follows. When 
we make our interactions longer-ranged in real space, we 
are, at the same time, making them shorter-ranged in 
momentum space. Let us start with a density-density 
interaction of the general form, Ylrr' ~ r')n(r)n(r'), 
where n(r) = J2 a c tr(. r ) c <?( r )- It turns out that the initial 
couplings in the effective low-energy field theory depend 
on the Fourier components V(q) of this interaction at 
q = and q = 2K; we will label these components as Vo 
and V2K, respectively. In the case of a long-range interac- 
tion V » V 2 k, and said interaction will therefore mostly 
cause forward scattering of the electrons while, in the case 
of a short-range interaction, V2K will become comparable 
to Vo , and thus back scattering begins to become compa- 
rable to the forward scattering. We would therefore say 
that we obtain a nematic phase when forward scattering 
dominates, while we obtain an antiferromagnetic phase 
when back scattering becomes comparable to the forward 
scattering. 

Interestingly, if we make the on-site (repulsive) inter- 
action weaker than the nearest-neighbor (repulsive) in- 
teraction, we may also obtain a quantum spin Hall phase 
for intermediate ranges of the interaction; whether or not 
it appears depends on the sign of one of our coupling con- 
stants, gE K , which we define in the main text, for a g iven 
value of £. This may be seen from Fig. 5 in Ref. [27t We 
were, however, unable to obtain a quantum anomalous 
Hall state in this case with any density-density interac- 
tion. 

Motivated by experimental datai^ and by our above 
conclusions, we then turn our attention to investigating 
the effects of an applied magnetic field on the antiferro- 
magnetic phase. A similar calculation has already been 
performed^, in which a self-consistent mean field ap- 
proach is employed, and a second order parameter corre- 
sponding to a "staggered spin current" order— is explic- 
itly introduced. Our variational mean-field analysis only 
explicitly includes the antiferromagnetic order parame- 
ter. We do this because, while at B = 0, the two order 
parameters are distinct in that they have opposite par- 
ity under mirror reflections, a finite external magnetic 
field, which is an axial vector, breaks the mirror sym- 
metry. Therefore, at B ^ 0, the two order parameters 
automatically mix; they belong to the same {A\ u ) repre- 
sentation of ^e, the point group for B ^ 0. Indeed, our 
variational wave function, which includes only the AF 
order parameter explicitly, leads to a finite expectation 
value of the "staggere d sp in current" order parameter 
when B =/= 0. In Ref. [28|, the development of a finite 



expectation value of the "staggered spin current" oper- 
ator was attributed to "the emergence of the n = 0, 1 
Landau levels (LLs) and the peculiar property of their 
wave-functions to reside on only one sublattice in each 
valley" . Here, we extend the argument and show that 
it must be present on more general grounds, and is not 
tied to the properties of the Landau levels. Note that, 
in this case, we are using mean-field methods to deter- 
mine the phenomenology of the broken-symmetry state 
in a case where we already know which symmetries are 
broken, based on our RG calculations, rather than as a 
means of predicting the broken-symmetry state. In fact, 
since the AF state is gapped, we expect that an expansion 
around the mean-field state should have a finite radius of 
convergence. Another reason why we employ mean-field 
methods rather than RG in this calculation is because, in 
the presence of a perpendicular magnetic field, the non- 
interacting energy spectrum is discrete and momentum k 
is not a good quantum number, thus making an analysis 
of the type used previously more difficult. 

In order to obtain a fit to the experimentally-measured 
gapii, we calculate the energy required to create a 
(charge neutral) particle-hole excitation, from which we 
can find the energy gap in the system. We find that it is 

m* 

E cx = mm[E ni + E n2 + —g*ui c (T 1 s 1 - r 2 s 2 )], (3) 

47T 

where tu c = eB/m*c is the cyclotron frequency for the 
effective mass m* , and 



E n = ^n(n-\)ul + k{Bf. (4) 

r and s are valley and spin indices; the former is ±1 for 
the ±K valley, and the latter is ±1 for spin up (down). 
Note that, if nj = or 1, then TjSj is locked to +1 for j — 
1 (particle) and to —1 for j = 2 (hole). This is simply the 
sum of the energies of the two states in the single-particle 
"auxiliary spectrum" , plus a non-universal term linear in 
the applied magnetic field, which depends on the coupling 
constant g* ~ V0—V2K discussed further in the text. The 
gap in our system is simply the minimum energy required 
to create such a particle-hole excitation. As discussed in 
the next paragraph, while the field dependence of the AF 
order parameter A(_B) is universally determined by A(0) 
through Eq. ([5]), the energy gap is not. The B-linear term 
has a coefficient that depends on the interaction in the 
combination, Vo — V2K, and thus the slope of the energy 
gap at high fields is independent of the zero-field value 
of the order parameter Ao = A(0). These results are in 
agreement with those found earlier in Ref. [28l 

In our method of determining A(_B), which was subse- 
quently reproduced in Ref. [28l we eliminate the coupling 
constants in the problem by rewriting the self-consistent 
equation in terms of Ao- This allows us to send the en- 
ergy cutoff in our problem to infinity. By doing so, we 
can obtain the scaling equation for the order parameter 
at a finite magnetic field whose dependence on coupling 



constants enters entirely through Aq: 



where 



F(a) = 1(a) + 



In 1 



(5) 
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and 1(a) is stated in Eq. (|C25[) . but may be very accu- 
rately approximated for any value of a as 



1(a) 



[( l/(0))-2/3 + 4.62/3 a 2]3/2' 



(7) 



and 1(0) w —0.0503 is the exact value of 1(a) at a = 0. 

While the above equation must, in general, be solved 
for A(B) numerically, it is possible to obtain approxi- 
mate analytic solutions in certain limits. For low fields, 
the order parameter has a quadratic dependence on the 
magnetic field, 



A(B) = A c 



8A ' 



(8) 



while, at high magnetic fields, the dependence on the 
field is 



A(B) = 



ln(w c /A ) + C" 



(9) 



where C ~ 0.67. This may be experimentally hard to 
distinguish from linear dependence. 

We find that E ex , given by Eq. ((3J, has a slight non- 
monotonic behavior for small fields. For example, if we 



choose 



to fit the experiment in Rcf. 0, we find a minimum of 
about 1.91A , occuring at around B = 0.047 T. At 
larger fields, we observe a "kink" , marking a transition 
to an approximately linear behavior, at a field of about 
0.45 T. This may be seen in Fig. [5j 

We also considered the high-field data for the v = gap 
given m Ref.fl In this case it is unclear from the low-field 
data what is the value of the energy gap at zero applied 
field, if any. Nevertheless, because of the weak logarith- 
mic dependence on Ao shown in Eq. ([9]). it is reasonable 
to fit this data using the expressions presented above and 
the same Ao = 0.95 meV as for the fit to the data of Rcf. 
4. In this case, we still obtain a non-monotonic depen- 
dence on the field, but with a very shallow minimum at 
about B = 0.017 T. In this case, we would need to as- 
sume Vq < V2K, which would require a non-monotonic 
density-density interaction. 

The rest of the paper is organized as follows. In Section 
II, we state the starting microscopic lattice Hamiltonian. 
In Section III, we derive the corresponding low-energy 
effective theory and the relations between the coupling 
constants in the low-energy theory and the interaction 
in the microscopic lattice Hamiltonian. Section IV is 
dedicated the determination of the phase of the system 



ff 



0.44 and A = 0.95 meV, which we use 




FIG. 1: (a) Illustration of the honeycomb bilayer lattice with 
only nearest-neighbor hopping terms. The circles represent 
the bottom layer, or layer 1, while the squares represent the 
top layer, or layer 2. The intralayer hopping is t, while the 
interlayer hopping is t±. We label the dimerized sites as a± 
and d2, while the non-dimerized sites are labeled 61 and &2- 
(b) Illustration of the (reciprocal) k space, with the parabolic 
degeneracy points K = g ^ x and K' = — K, where a is the 
(real-space) lattice constant, marked. 



as a function of the range of the interaction using the 
results of Section III as the initial conditions for the RG 
analysis detailed in Ref. [27]. In Section V, we present our 
variational mean field analysis of the antiferromagnetic 
state in the presence of an applied magnetic field and 
calculation of the energy gap within this approximation. 
We present our conclusions in Section VI, and provide 
mathematical details in the appendix. 



II. MICROSCOPIC LATTICE HAMILTONIAN 

Our starting point will be a tight-binding A-B stacked 
honeycomb bilayer lattice with only nearest-neighbor 
hopping terms and a two-particle interaction. The lattice 
and the corresponding reciprocal momentum (k) space 
are illustrated in Figure [T] By including only nearest- 
neighbor hopping, we are fine-tuning our system so that 
it only possesses two parabolic degeneracy points, rather 
than four Dirac cones. As stated earlier, this will make 
our system unstable to symmetry-breaking phases for ar- 
bitrarily weak interactions, thus justifying the use of pcr- 
turbativc RG methods. The Hamiltonian for this system 



H = + Hq + Hj + 
where the kinetic energy terms^ are 



(10) 



#0 ="*E E [bl(B.+S)a llT (R)+bl(R-S)a 2a (K)+h, 

R,<5 <r=T,4. 

(11) 

and 

H o =-«x£ £ [al(K)a 2a (K) + h.c.]. (12) 

R <t=t,4- 



Here, t and t± are the intra-layer and inter-layer nearest- 
neighbor hopping integrals, respectively, S is a vector that 
connects an a\ site to a nearest-neighbor b\ site. The 
position vectors R that we sum over are the positions of 
the dimerized sites. The three possible values of 8 are 
■^ax + iay, — ^ax + \ay, and —ay. Whenever there 
is a sum on <5, we sum over all three vectors, while, if S 
occurs without a sum over it, then we choose one such 
vector. 

The interaction terms are given by 



H i = \ E E V \\ ( r - r ')K(r) - - !] ( 13 ) 



k— 1 rr' 



and 



ff/ = J2 - r ')Mr) - l][n 2 (r') - 1]. 



(14) 



Here, r runs over all projections of the position vectors of 
the lattice sites onto the plane of the sample, and thus is 
entirely in the xy plane. Cfc CT (r) is the annihilation opera- 
tor for a particle at site r, and rifc(r) = c lcr( r ) c fco-( r )- 
The interaction V(r) is assumed to depend only on dis- 
tance, i.e., V(r) = V(|r|). For convenience, we intro- 
duced the notation Vji(r) = V(r) and Vj_(r) = V(r±cz), 
where c is the distance between the two layers. The sys- 
tem represented by this Hamiltonian will be at half filling 
when the chemical potential is zero. This follows from the 
fact that the Hamiltonian is invariant under the particle- 
hole transformation, 



«L(R)> 



oi ff (R) = 5 t 1[T (R) 

02<r(R) = 

6 lff (R) = -&L(R), 

6 2CT (R) = bl(R). 



(15) 
(16) 
(17) 
(18) 



Using this, the calculated expectation value of the parti- 
cle number on a given site can be shown to be 1. 



III. LOW-ENERGY EFFECTIVE THEORY 

Before performing our RG analysis of the above Hamil- 
tonian, we first derive the low-energy effective theory as- 
sociated with it. We assume that the kinetic energy terms 
are dominant, and that the interactions are small. This 
will allow us to focus on the interaction-induced scat- 
tering processes in the vicinity of ±K = ±47r/(3v / 3ct)x. 
Note that there are four sites per unit cell, and therefore 
there are four bands. For every state with energy E(k), 
there is another with energy —E(k). To proceed, we first 
write the partition function for the system as a coherent- 
state path integral and integrate out the dimerized sites 
ai and a.2- To see that the high- energy modes are as- 
sociated with these sites, note that the separation of the 
high-energy, split off, bands is set by the hopping integral 



t± between a\ and a 2 . We then expand around the two 
parabolic degeneracy points in the Brillouin zone, asso- 
ciated with the remaining two bands. The derivation of 
this theory is a straightforward, if tedious, generalization 
of the steps followed in Ref. 0. The low-energy effective 
theory that we obtain is 



Z = J £>[V>*,V>] exp (-J^drL^j , 
where the Lagrangian L c ff is given by 



(19) 



1 + 







L cS = Jd 2 K^ 

+ 1^2 3s J d 2 R(^SiP) 2 -i£ I d 2 R^^. 



(20) 

Here, V( r > T ) = [V>t( r > T )> i>±( T > T )F is an ei § nt_ 
component spinor in layer (1, 2), valley (K, — K), and 
spin (f, I) space, and 



i/'a 



' 6i,K,ff(r,r) ' 
&2,K,ff(r,r) 

bl,-K,a( r > T ) 
&2,-K,cr(r,T), 



(21) 



The matrix, H(p), is given by 



H(P) 



Px~P 



2m* 



*E 3 



PxPy , 

m* 



(22) 



where m* = 2t±/9a 2 t 2 is the effective mass, Y, x = 
l 2 a x l 2 — 72 ® I2, and E a = r z a y l 2 = 71 ® h- In 
the former definitions, the t matrices act in valley space, 
the a matrices act in layer space, and the third matrix, 
which is the identity in both cases, acts in spin space. 
The sum on S is over all 16 matrices of the form, Ti<jjl 2 . 
These 16 matrices may be classified according to the rep- 
resentation of the space group under which our system is 
symmetric, namely the D^d point group plus appropri- 
ate translations, that they transform under. A complete 
classification of the 16 possible 4x4 matrices that act in 
valley and layer space is done in Ref. H3; we repeat those 
results here for convenience: 
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A 1K + 
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T y <J X 


A 2K - 
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T y a y 


E K + 


■ (r x h 


-r y a 



The ± at the end of each representation denotes how the 
associated matrices transform under time reversal. The 



coupling constant gs is that associated with the repre- 
sentation that the matrix S belongs to. Since the sum 
includes all 16 matrices, there are nine coupling constants 
in all. It can be shown^ that there are only nine inde- 
pendent coupling constants due to the symmetries of the 
underlying lattice and Fierz identities. In the notation of 

r> c □ ( c ) ( c ) ( c ) ( c ) 

Ret. |6J, 9A lg = g Al , gA 2g = 9d 2 ' 9e 9 = g A2 , 9a 1u = g B ^ 



9 A 



o (c) 



9A 2 



<h 



c) 



and 



(c) fc) 

: ffci, 9e u = 9b[, gA 1K 

(c) 

9e k =9p- 

Provided that we start with density-density interac 
tions only, as we do in this case, the only (initial) non 
zero coupling constants in this theory are 



where 



9A lg 
9a 2u 
9e k 



v ±iN 

V\\,2K 



K^ll.o 



V±, N )A 
V, N )A 

U C •■ 



■gV\\,2K Aic, 
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R,<5 

J2 V \\( R ) cos(2K-R), 



(23) 
(24) 
(25) 



(26) 
(27) 
(28) 
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We note that gA lg and gA 2u depend only on q = Fourier 
components of the interaction, and thus we may say that 
they give us a measure of the strength of the forward 
scattering induced by said interaction. Likewise, we note 
that gE K only depends on the q = ±2K Fourier com- 
ponents, and thus it gives us a measure of the strength 
of the back scattering. Furthermore, we see that £U 2 „ 
depends on the difference between an intra-layer interac- 
tion and an inter-layer interaction, and thus it may be 
seen as a measure of the imbalance between these two 
interactions. 

Note that our theory includes a quadratic, chemical 
potential-like, term. We may think of the undetermined 
constant fjf as being chosen in such a way as to can- 
cel out the quadratic terms that are generated from the 
quartic terms under RG. We require that this occur be- 
cause we know that our original lattice model is at half 
filling (that is, it possesses particle- hole symmetry), and 
therefore this must be reflected in our effective low-energy 
theory as well. 

We could, in principle, have also determined the value 
of // when we wrote down the above effective low-energy 
theory. Strictly speaking, we should not simply drop 
all of the modes above the cutoff, as we did here, but 
rather integrate them out in a pcrturbativc scheme sim- 
ilar to what is done in an RG analysis. This, at first 
order, will not change our quartic terms because it only 
generates the tree-level quartic terms. However, it will 
generate both tree-level and one-loop contributions to 
the quadratic terms. It would, however, be somewhat 



cumbersome and, given the above particle-hole symme- 
try arguments, equally unnecessary, to determine these 
one-loop contributions to the chemical potential. 



IV. DETERMINATION OF PHASES 

We are now ready to describe the results of our RG 
analysis. We consider two forms of the microscopic in- 
teraction, which are given by Eqs. (JTJ) and ([2]). In both 
cases, we determine our initial couplings by first using 
Eqs. (|23|) -([25 |) to determine the ratios, gA 2u /9A lg and 
9e k I 9A lg i as a function of the range of the interaction. 
We then consider different values of g Al to multiply 
these ratios by; we may therefore think of this value of 
gA lg as determining the overall strength of the interac- 
tion. We then use these couplings as the initial con- 
ditions for the RG equations derived in Ref. [13, which 
we integrate numerically. We then adjust the temper- 
ature until we see the couplings diverge as we take the 
scale parameter to infinity, i.e., until we reach the critical 
temperature. We monitor the susceptibilities to differ- 
ent phases as we lower the temperature, and consider a 
phase to be present below the critical temperature if the 
corresponding susceptibility diverges as we approach the 
critical temperature. 



A. Screened Coulomb-like interaction; two-plate 
case 

The first interaction form that we will consider is a 
screened Coulomb-like interaction, with the screening be- 
ing due to two infinite planar conducting plates, between 
which the charge is located. We consider this case be- 
cause of experiments with gated bilayer graphene, such 
as the experiments in Refs. 13 and 0; in these cases, the 
gates will serve as the conducting plates. We will assume 
that the distance between the plates is much larger than 
the distance between the two layers so that we may as- 
sume that the particles are exactly halfway between the 
two plates. If the distance between the plates is £, then 
the interaction is given by 



V(r) = U J2 



(-1)" 



(29) 



As shown in Appendix [Al for |r| S> £, we may approxi- 
mate the sum as 



V(r) « U, 



2v / 2e~ 7r l r l^ 



(30) 



The above form is useful in practice when evaluating the 
values of the initial coupling constants. In the opposite 
limit, |r| -C £, we may simply approximate the interac- 
tion with the first few terms of the sum around n = 0. 

Note that, as is, the on-site interaction given by our for- 
mula is infinite, and thus it would give us infinite values 
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FIG. 2: Phase diagram for the two-plate screened Coulomb- 
like interaction, Eq. (|29[) . as a function of the dimensionless 
coupling strength ! j^-gA lg and the interaction range £ in units 
of the lattice spacing a. For very short ranges, we find an 
antiferromagnetic (AF) phase. As we increase the range, we 
enter a region where the susceptibilities toward both the AF 
and nematic (N) phases diverge as we lower the temperature. 
Determining whether or not these two phases truly coexist 
requires a theory valid below the critical temperature, the 
development of which is beyond the scope of the present work. 
Further increasing the range brings the system into a pure 
nematic state. Note that the critical range for each of these 
transitions is weakly dependent on the coupling strength, and 
corresponds to effectively turning off back scattering. 



FIG. 3: Phase diagram for the one-plate screened Coulomb- 
like interaction, Eq. (|31|l . as a function of the coupling 
strength <?a 1s and the interaction range £. The phase dia- 
gram is qualitatively identical to that obtained for the two- 
plate case, except that the critical ranges are smaller. 



happens when £ exceeds about 10 lattice spacings. Note 
that there is a weak dependence of these critical ranges 
on the initial value of gA ± ■ 



B. Screened Coulomb-like interaction; one-plate 
case (dipole-like interaction) 



for the initial coupling constants. However, we recognize 
that, for two electrons on the lattice that are sufficiently 
close, the electrons are not localized at a single point, 
but rather their wave functions have a finite extent in 
space. This will render the on-site interaction finite. As 
a simple model of this effect, we set the on-site interac- 
tion equal to some constant A times the nearest-neighbor 
interaction. In our calculations, we set A = 1.2. 

The resulting phase diagram is shown in Figure [2j We 
see that, for very short ranges of the interaction £, the 
system is in an antiferromagnetic state, in agreement 
with the zero-temperature results obtained in the pre- 
vious work*. As we increase the range, we will enter a 
region in which both the AF and nematic susceptibilities 
diverge^!. This happens when £ is larger than about two 
lattice spacings. This indicates a possible coexistence of 
the two phases. However, to determine if such a coex- 
istence is in fact present would require a theory valid 
below the critical temperature. The construction of such 
a theory is beyond the scope of the present work. If we 
increase the range even further, we enter a pure nematic 
state, again in agreement with the previous work^. This 



The second form of the interaction that we consider is 
a dipolc-likc interaction much like the one produced by 
an electron in the presence of a single infinite conducting 
plate. This interaction has the form, 



V(r) = U 



1 



1 



(31) 



This interaction has a longer range than in the two-plate 
case, since this falls off as r~ 3 for long distances, rather 
than as an exponential. As in the previous case, this 
formula, as is, will give us an infinite on-site interaction. 
We use the same method as before to render this inter- 
action finite, and we again set A = 1.2. The resulting 
phase diagram is shown in Figure [3] We note that it is 
qualitatively identical to that obtained from the previous 
case, except that the critical ranges are smaller. We enter 
the AF + nematic "coexistence" region when £ exceeds 
a value between 0.4 and 0.6 lattice spacings and the pure 
nematic state when £ exceeds a value between 4 and 6 
lattice spacings, depending on the initial gA lg ■ 

Note that, throughout this section, we have been 
working with monotonically-decreasing repulsive density- 



density interactions, and thus we only observe two of the 
possible phases that we may find in the system. For 
short-range interactions, we start in the upper-right-hand 
corner of Fig. 5 of Ref. [27|, corresponds to the AF phase. 
As we increase the range, we move toward the center of 
the diagram, while staying within the upper-right quad- 
rant, passing through the AF + nematic "coexistence" 
region and ending in the pure nematic region. While, 
for A = 1.2 we only find the AF and nematic phases for 
the one- and two- plate cases that we considered, there 
are more possible phases, even for density- density inter- 
actions. For example, for A < 1, such that the on-site 
repulsion is weaker than the nearest-neighbor repulsion, 
and thus the repulsive interaction is non-monotonic in 
real space, the initial value of gE K may become negative. 
Under such conditions, we will find that the suscepti- 
bility toward the quantum spin Hall phase will diverge, 
though never along with the antifcrromagnetic suscepti- 
bility. This is illustrated in Fig. 5 of 



V. ANALYSIS OF THE 
ANTIFERROMAGNETIC STATE IN AN 
APPLIED MAGNETIC FIELD 

A. Field dependence of the AF order parameter 

We now turn our attention to the effects of a magnetic 
field on the antifcrromagnetic state of the system at zero 
temperature. This investigation is motivated by the fact 
that a gap is observed in some experiments^ and by 
the fact that we predict such a phase for short-range in- 
teractions within RG. We will investigate these effects 
within the framework of variational mean field theory. 
We employ this method, rather than RG, because, in 
the presence of a perpendicular magnetic field, the non- 
interacting energy spectrum for our problem is discrete, 
rather than continuous, and thus an RG analysis of the 
type employed above would be more difficult. We start 
by writing down a Hamiltonian corresponding to our ef- 
fective low-energy field theory and introducing the orbital 
effects of the magnetic field via minimal substitution. We 
will neglect the Zceman effect in this case, since the spin 
splitting (m* /m e )g/iBB ~ (m* /m e )u> c is small compared 
to the orbital splitting uj c . This Hamiltonian is 

H = [ d 2 r^{r) I J2 d a (Tv)T, a J V(r) 

+ IE / d 2 rg s [^(r)Sm} 2 , (32) 
s J 

where tv = p — -A, A is the applied magnetic vector 
potential, 

*M-^, (33) 

and the sum on S stands for the four-fcrmion interac- 
tion terms that appear in Equation (|2"U1) . Here, ip is a 



field operator, not a Grassman field as it was in previous 
sections. 

Our variational mean field calculation will proceed as 
follows. We start by adding and subtracting a source 
term for the antiferromagnetic order parameter, 



A / d 2 r^(r)l 2 a z s z i;(r), 



(34) 



in the Hamiltonian. We now define two parts to the 
Hamiltonian, a "non-interacting" part, Hq, 

Ho = J dW(r) ( £ d°(7r)S°) tf(r)+ 

\a=x,y / 



+ A / d 2 r^(r)l 2 a z s z ^(r), 



(35) 



and an "interaction", Hi, 

H i = IE / d 2 rg s [^(r)Si>(r)F- 

S,T 

- A J d 2 r^{Y)l 2 (j z s z ip{Y). (36) 

We then exactly diagonalize Hq, find the expectation 
value of the full Hamiltonian with respect to the ground 
state of Hq, and minimize the result with respect to A. 
We will provide the details of the diagonalization in Ap- 
pendix [B] and simply quote the main result here. In the 
Landau gauge, our states can be described in terms of 
a wave number k, an orbital index n, a valley index r, 
and a spin index s; the corresponding energy eigenvalues 
are E n = ±y/n(n — T)w| + A 2 , where ui c = -^f- is the 
cyclotron frequency of the electrons. For n > 2, each of 
these energy levels is four-fold degenerate due to valley 
and spin degeneracies. The levels given by n = and l 
are also four-fold degenerate, but this time due to orbital 
(n = or l) and spin degeneracies. 

We may now rewrite the field operators in terms of 
these eigenstates: 

k,n,r,s 

(37) 

where i/j^n I s represents the positive (negative) energy 
state for a given wave number and orbital, valley, and 
spin indices. We have chosen our operators a and b such 
that they annihilate the ground state |0). We may now 
take the expectation value of our Hamiltonian with re- 
spect to the ground state of H and minimize the result 
with respect to A; we will provide some details on how 
this was done in Appendix [Cj Upon doing so, we find 
that the equation for A is, assuming A > 0, 



> N 

E 



where 



9cS 



4tt 



(9A lg +9A 2u +4g EK ) 



(38) 



(39) 



and N is an upper cutoff on the orbital index, which we 
impose because our theory only works for low energies 
and because, in reality, we do not have electronic states 
in our system at arbitrarily large energies. In the limit 
of zero magnetic field, we may treat the sum as a Rie- 
mann sum, with e = nuj c and 5e = uj c ; our equation then 
reduces to 



de 



A 2 



= 1. 



(40) 



where f2 is an upper cutoff on the energy and we intro- 
duce Ao as the value of the AF order parameter in the 
absence of an applied mag netic field; N is related to the 
energy cutoff by ft = lj c WN(N — 1). We have verified 
this result with a separate calculation. As is, we cannot 
send the upper cutoff to infinity in our equation without 
encountering a divergence in the sum, since the summand 
only decreases as —. We can, however, rewrite the equa- 
tion in such a way that we can do this if we eliminate g c ff 
in favor of Ao- We give the details on how we do so in 
Appendix E2 the final result is 



1(a) + In 

a 



1 



« 2 + ! 



(41) 



where a = A(B)/co c , (3 = w c /A , A(B) is the order pa- 
rameter at finite field, and lu c is the cyclotron frequency 
of the electrons in our system. The function, 1(a), the 
exact form of which we state in Appendix [C] is monoton- 
ically decreasing, falls off as a~ 3 for large a, and may be 
very accurately approximated for any value of a by 



1(a) « - 



[(_lj(0))-2/3 + 4.62/3 a 2]3/2^ 



(42) 



where 1(0) w —0.0503 is the exact value of 1(a) at a = 0. 

While the above equation must, in general, be solved 
numerically, we may obtain analytic expressions in two 
limiting cases, namely the large and small j3 limits (cquiv- 
alently, for large and small applied magnetic fields). Let 
us first consider the large (3 limit. This means that the 
right-hand side of our equation is large and positive, and 
thus, as implied by our above discussion, a should be 
small. In this limit, we may set the first and third terms 
on the left to their values at a = 0, since both are fi- 
nite at this point, while the second term diverges. Our 
equation becomes 



1 



C = ]n/3, 



(43) 



where C is the value of the first and third terms at a = 0; 
its value is approximately 0.67. Here, we are assuming 
that a is positive. Solving for a, we find that 



1 



a 



ln/3 + C" 



or 



A 



ln( Wc /A ) + C- 



(44) 



We see that the behavior is almost linear in the magnetic 
field, but with a logarithmic correction. 

Next, we consider the small /3 limit. In this case, the 
right-hand side of Eq. (|3"5)) becomes large and negative, 
and thus a should become large. We would be tempted 
to drop all but the third term, since the first two terms 
go to zero for large values of a. This would give us a 
result that is only accurate at constant order in /3, how- 
ever. This is because, if we expand the third term in 
powers of a -1 , then the next-lowest order term after the 
logarithmic term is of order a -1 , and the second term in 
our equation is also of this order; in fact, it cancels this 
term exactly. We may still drop the first term, since, as 
stated above, the lowest-order term that it contributes is 
of order a -3 . In this case, our equation becomes 



- - In ( 1 + Jo 2 + 3 

a 



= In/?. 



(45) 



Again, we are assuming that a > 0. If we take the expo- 
nential of both sides, we get 



e' 1 ^ { l + A /a* + 



2 i 3 

4 



1 



We will now expand the left-hand side in powers of a^ 1 
to the order a~ 2 . Doing so, but first pulling out a factor 
of a from the second factor, we get 



a 1- 



1 

8^ 



1 



If we rearrange this, we finally arrive at the quadratic 
equation, 

8/3a 2 -Sa-P = Q. 

If we solve this equation and take the positive solution, 
we get 



1 



1 



\fi 2 



2/3 



1 





Rewriting this result in terms of A and w c , we get 



A = A 



8A ' 



(46) 



We see that, for low fields, the antiferromagnetic order 
parameter increases quadratically with the field. 

We now solve Equation (|41[) numerically; the numer- 
ical result, along with the low- and high- field limits de- 
rived above, is plotted in Figure |4j If we look at our 
low- and high-field expressions, we sec that the slope of 
our low-field approximation increases with B 7 while our 
high-field approximation has a decreasing slope. This 
implies that there should be a maximum slope to the 
exact curve. We determined the maximum slope of the 
A(_B)/A versus oj c / A curve, and found that it is about 
0.2681, and occurs when ui c / Aq ~ 2.432. These values 
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FIG. 4: Plot of the solution to Equation (jiTj) . The solid line is 
the numerical solution, while the dashed lines are the solutions 
in the low- and high-field limits. The vertical axis is the value 
of the order parameter A divided by its zero- field value Ao. 
The bottom horizontal axis is the cyclotron frequency ui c = 
eB/m*c divided by Ao, while the top horizontal axis is the 
applied magnetic field B divided by Ao; in determining the 
latter from the former, we assumed that the effective mass^ 
m* = 0.028m e , anc£ A = 0.95 meV. 



are independent of the values of m* and Ao- Using the 
experimentally-determined valued of the effective mass, 
to* = 0.028m e , and the experimentally-determined value 
of the order parameter at zero field^, Ao = 0.95 meV, we 
may determine the maximum slope of the A(i?) versus 
B curve. We find that the slope is 1.11 and that it 
occurs at a field of about 0.56 T. 



B. Comparison with experimental results 

We would now like to compare our theoretical results 
to the experimental datai. First of all, we note that A 
is not the energy gap in our system. In fact, the energy 
eigenvalues stated earlier are an "auxiliary spectrum" , 
and do not represent the true (many-body) energy spec- 
trum of our system. As an approximation to the actual 
energy gap, we will consider particle-hole excitations of 
the "vacuum" , or trial ground state, for our system. We 
construct a state, a^bp |0), where a and f3 stand for the 
full sets of quantum numbers describing the particle and 
hole states, respectively, and find the difference between 
the expectation value of our Hamiltonian for this state 
and that for the trial ground state. For both states, we 
assume the value of A that is obtained from the mini- 
mization condition, Eq. (|4ip. The states a and /3 that 
result in the lowest value of the excitation energy will be 
taken to give the actual energy gap. We will provide de- 
tails of the derivation in Appendix [D] and simply quote 
the final result here. If we take the electron (hole) state 
a (/3) to have an orbital quantum number ri\ (712), wave 



number q\ ((72), valley index n (T2), and spin S\ {S2), 
then the excitation energy is 

TO* 

E cx = E ni + E n2 + —g*u! c (T 1 s 1 - t 2 s 2 ), (47) 

47T 

where g* = gA lg + 9A 2u ~ ^9e k - Here, r k = ±1 if the 
state exists in the ±K valley, and sj, = +1 ( — 1) for 
a spin up (down) state. The single-particle energies are 
E n = \Jn{n — 1)cj2 + A 2 . Note that this energy does not 
depend on the wave numbers of the states, and only de- 
pends on the valley and spin indices via their products. 
If one of the states is one of the lowest Landau levels, 
given by n k — or 1, then these products are locked to a 
specific value. To be exact, if n\ = or 1, then t\Si = 1. 
Similarly, if n% = or 1, then T2S2 = —1- We see that 
this energy includes a term linear in the magnetic field, 
and is in agreement with the results obtained in Ref. l28l . 
The key difference between our derivation and that pre- 
sented in Ref. is that we did not need to assume the 
presence of another "order parameter" (in fact, as we will 
explain shortly, this other paramter is not really an or- 
der parameter in the sense that it breaks any additional 
symmetries), corresponding to the matrix t z \is z in the 
notation of the present paper, or a "staggered spin cur- 
rent" stated to obtain this linear term, assuming that 
we properly calculate the excitation energy (i.e., we cal- 
culate it from the full Hamiltonian rather than assume 
that the gap in the single-particle "auxiliary spectrum" 
is the observed gap). In fact, the above result would 
not have changed had we had included this parameter in 
our variational analysis, assuming that it was sufficiently 
small — it would have only introduced a constant shift 
to the energies in the "auxiliary spectrum" and left the 
associated wave functions unchanged. 

We present a plot of part of this excitation spectrum 
in Fig. [5l We find that the gap for low fields is, in fact, 
not given by taking n\ and n 2 to both be either or 
1 , which would correspond to the lowest-energy states in 
the "auxiliary spectrum" . Instead, it is given by taking 
ni = ri2 = 2, with TiSi = —1 and T2S2 = 1. At higher 
fields, however, we find that excitations with ri\ and n 2 
both equal to cither or 1 do, in fact, give us the actual 
gap. To obtain the value for g* , we fit the slope of our pre- 
dicted high-field gap at around 2.5 T to the slope found 
in Ref. of 5.5 Assuming that the effective mass is 
given by the experimental valued of to* = 0.028m e , where 
to c is the mass of an electron, we obtain ^g* = 0.44. 
Note that this differs slightly from the value used in Ref. 
[28l . namely ^g* = 0.4; this is the value that we would 
obtain if we instead fit the slope of the high-field gap at 
the point where the AF order parameter reaches its max- 
imum slope to the experimental value. For the value of 
g c s that we use, we find that the gap has a minimum at 
a non-zero value of the field; the minimum is reached at 
a field ofB»i 0.047 T, and is E CJ w 1.91A . We also find 
a kink in the field dependence of the gap, at which the 
gap goes from being given by ni = n 2 = 2, TxSi = — 1, 
and T2S2 — 1 to being given by n\ and 712 either or 1. 



This kink appears at a field of B 0.45 T. 
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be observed in the experiments due, for example, to the 
fact that, at finite temperature, any sharp features that 
would have appeared at zero temperature are "washed 
out" , thus introducing uncertainty into any energy gaps 
extracted from the data. It is also possible that these 
features, which are predicted from a mean- field calcula- 
tion, will be removed once fluctuations are taken into ac- 
count. The development of a more sophisticated method 
for treating this problem is therefore of interest, but it is 
beyond the scope of the present paper. 



Symmetry analysis in the presence of an 
applied magnetic field 
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FIG. 5: Plot of the excitation spectrum, Eq. (|4"T)>. The verti- 
cal axis is the excitation energy in units of Ao, and the hor- 
izontal axes are the same as in Fig. [4] Each curve is labeled 
according to which electron and hole states are occupied, with 
a label of the form (electron, hole), with the orbital index n 
and the sign of the product of the valley and spin indices, rs, 
indicated. 

We also considered the data from Ref. [l|. In this case, 
we fit our expression to the v = gap presented therein. 
Note that, from the low-field data, it is unclear what the 
size of the energy gap, if any, is. Nevertheless, we can 
still obtain a fit to the slope of the gap at high fields, 
since the zero-field value of the order parameter only en- 
ters via a small logarithmic correction to the high-field 
slope. The experimentally-determined slope is 1.7 
If we perform our fit in the same way as before, we ob- 
tain ^-o* = —0.018. In this case, since we obtain a 
negative value for g* , we will find that energy of the n\ 
and n 2 equal to or 1 excitation has a non-monotonic 
dependence on the magnetic field. In fact, it gives us the 
energy gap for all fields. It possesses a very shallow mini- 
mum of 1.99985A at a field of about 0.017 T and, unlike 
the previous case, there is no kink. In this case, we can- 
not completely rule out the possibility of the gap actually 
possessing such a minimum on the basis of the data given 
in Ref. [l] alone, due to the lack of data at low fields. Note 
that we required a negative value of g* , which would im- 
ply that AgE K > gA lg + 9A 2u i to fit the data. Satisfying 
this inequality would require either an attractive inter- 
action or one that is non-monotonic; this may be seen 
by noting that it is equivalent to V\\ 2 k > M|,0i which 
cannot be satisfied for any monotonically-dccreasing re- 
pulsive interaction. 

Note that, while we are able to fit the experimental 
data^ at high fields, we also predict finer features at low 
fields that are not resolved in these experiments, namely 
a slight non-monotonic behavior of the gap and, in the 
case of our fit to the data from Ref. 0, a "kink" . It is pos- 
sible that such features are, in fact, present, but cannot 



In the absence of a magnetic field, the honeycomb bi- 
layer lattice considered in this paper possesses a 
point group symmetry. The AF and "staggered spin cur- 
rent" orders transform under different representations of 
this group — the AF state transforms under the A 2u 
representation, while the "staggered spin current" state 
transforms under the Ai u representational. Physically, 
this means that the AF order parameter is even under 
mirror reflections and odd under C' 2 rotations, while the 
"staggered spin current" order parameter is odd under 
mirror reflections and even under C' 2 rotations. Note that 
both orders are odd under inversion. In fact, the two also 
transform differently under time reversal; the AF order 
parameter is odd, while the "staggered spin current" is 
ever*2£. This means that the expectation value of the 
"staggered spin current" operator must vanish in the AF 
state. 

When we apply a magnetic field, however, the point 
group is reduced to Sq . This is because the magnetic field 
is an axial vector that is odd under mirror reflections and 
even under inversion. In this case, the AF and "staggered 
spin current" orders transform under the same represen- 
tation, namely the A\ u representation^. Physically, this 
is because the mirror reflection and C' 2 symmetries are 
no longer present. As pointed out above, the two or- 
der parameters transform differently under time reversal; 
however, time reversal symmetry is broken in the pres- 
ence of a magnetic field. This means that the two orders 
no longer break different symmetries, and thus there is 
nothing preventing the system from acquiring a non-zero 
expectation value of one of these order parameters in the 
presence of the other. 

The development of a finite expectation value of the 
"staggered s pin current" operator was correctly pointed 
out in Ref. I28L but was attributed to "the emergence 
of the n = 0, 1 Landau levels (LLs) and the peculiar 
property of their wave-functions to reside on only one 
sublattice in each valley" . Here, we show that it must be 
present on much more general grounds, and is not tied 
to the properties of the Landau levels. 

At B = 0, the AF order parameter breaks time reversal 
and inversion symmetry, but it preserves mirror reflection 
symmetry^. Therefore, the wave functions for this state 



are eigenstates of the reflection operators as well, and 
may be classified as even or odd under them. Let us 
now consider the expectation value of an observable O 
that transforms under the A\ u representation of the D$d 
point group, such as the "staggered spin current" order 
parameter. This operator will have the property that 
any mirror reflection ad will anticommute with it, i.e. 
UdO = —Odd- Because of this, the expectation value of 
the "staggered spin current" operator with respect to the 
AF state of the Hamiltonian must be zero. 

Let us now consider the case in which B ^ 0. As 
stated before, this will break the mirror reflection sym- 
metry of our system. However, it is symmetric under 
such a reflection followed by a reversal of the magnetic 
field (i.e., B — > —B). This means that we may classify 
all eigenstates as even or odd under this combination of 
operations. This means that, in terms of the eigenstates 
of the Hamiltonian at B — 0, we may write the new AF 
state of the Hamiltonian in the presence of an applied 
field as 

\AF{B)) =Y,Wf\B)\i,e) + a\°\B)\i,o)], (48) 

i 

where \i,e) (\i,o)) is a general even (odd) state of the 
zero-field Hamiltonian. One set of the a coefficients (i.e., 

oif or otf 1 ) must be even functions of B, while the other 
must be odd. If we now calculate the expectation value 
of O with respect to this state, only matrix elements 
that mix states of opposite parity under reflections will 
appear: 



» 



\0\j,o) 



c.c. 



(49) 



Since one of either the af^ or must be even func- 
tions of B, while the others must be odd, we see that the 
expectation value of O must be an odd function of B. 

If we calculate the expectation value of the "staggered 
spin current" operator for the trial ground state that we 
work with above, we find that it is a linear function of B. 
This is consistent with our general conclusions and with 
the observation made by Kharitono\«2&. 



D. Corrections to the "auxiliary spectrum" 

There is one other point that we wish to address. 
Throughout this calculation, we have been assuming that 
the presence of the order parameter opens up a simple 
gap in the spectrum, modifying the low-energy disper- 
sion to have the form, E(k.) = ±^/[e(k)] 2 + A 2 , where 
e(k) is the dispersion in the absence of the order param- 
eter. This assumption is, in fact, not entirely true; in 
reality, the dispersion acquires a "Mexican hat" shape. 
This can be seen by taking the continuum limit of the 
effective action given by Eqs. (19)-(21) in Ref. 0, which 
will give rise to an additional term (vF<]/t±) 2 d/dT, where 
vp = |at, in Eq. (24) in the same work. In the absence 



of an applied magnetic field, but in the presence of the 
antifcrromagnetic order parameter Ao, the dispersion is 



E(k) = ± 



A 2 



i + (v F k/t ± y 



(50) 



If we expand this expression to fourth order in fc, we 
obtain 



E(k) w ±A 



k' 2 



2m* A 



(51) 

We see that the quadratic term is negative, while the 
quartic term is positive, thus giving this dispersion the 
"Mexican hat" shape. Note that it is the presence of a 
non-zero Ao that changes the scaling dimensions of our 
operators and necessitates the inclusion of the above- 
mentioned term for a complete description of the low- 
energy behavior. For Ao = 0, which is the starting point 
of our RG analysis, the scaling dimensions of our opera- 
tors are the same as they were in Ref. @, which justifies 
neglecting this term in our RG analysis. Performing the 
standard minimal coupling to an external vector poten- 
tial, we have looked at what effect this would have on 
the lowest Landau levels, and we find that it lifts the or- 
bital degeneracy of the n = and 1 Landau levels. To 
be exact, these two levels take the form, 



E n = 



1 + fMn- 



(52) 



For small magnetic fields, for which uj c <C we find 
that the splitting between these levels is approximately 
^A = aB, where a m 0.0106 Sf^. This is small com- 
pared even to the Zeeman splitting, which is also linear 
in B, but with a slope of about 0.232 Therefore, we 
may safely neglect this effect in our calculations. 



VI. CONCLUSION 

We have employed weak-coupling pcrturbative RG 
methods to investigate the different phases that fermions 
on a honeycomb bilayer lattice with finite-range interac- 
tions will enter. The use of these methods is justified 
since we only include nearest-neighbor hopping terms, 
resulting in a band structure with two quadratic degen- 
eracy points and therefore in a finite density of states at 
the Fermi level at half filling. We considered two forms 
of the interaction, a screened Coulomb-like interaction 
much like the one produced by a point charge situated 
exactly halfway between two infinite parallel conducting 
plates, as well as that produced by a point charge in the 
presence of a single conducting plate. For all cases, we 
determined what phase the system enters as a function 
of the range of the interaction. 

We found that the system, for both forms of the inter- 
action, enters an antifcrromagnetic state for short ranges 



and a nematic state for long ranges, in agreement with 
the previous work^. For intermediate ranges, we find 
that the susceptibilities towards both the antifcrromag- 
netic and nematic phases diverge, though not necessarily 
with the same exponent^. This indicates a possible co- 
existence of the two phases. To determine whether the 
phases truly coexist, or if only one appears, would require 
a theory that is valid below the critical temperature. The 
development of such a theory is a problem of great inter- 
est, but is beyond the scope of the present work. 

While we find that short-ranged interactions result in 
an antiferromagnetic state, which is gapped, while long- 
ranged interactions result in a nematic state, which is 
gaplcss, the ranges at which we see a transition from one 
to another are too short to explain why the experiment 
in Ref. observes an energy gap in their sample, while 
that of Ref. [U finds evidence for a nematic state in their 
sample. We note, however, that the sample studied in 
Ref. HI has a higher mobility than that studied in Ref. 0, 
which suggests that the degree of disorder present in the 
system is a major influence on the phase that appears. 
Throughout our work, we have assumed a clean sam- 
ple; the incorporation of disorder into our system may 
therefore be important in better explaining this appar- 
ent discrepancy in the experimental results. 

Another issue is the fact that we predict an antiferro- 
magnetic phase, which breaks a continuous SU (2) spin 
symmetry, which is impossible in two dimensions at finite 
temperature. We should therefore view the divergence of 
the susceptibility toward this phase as simply identifying 
the dominant ordering tendency in this case, even if there 
is no actual symmetry breaking, as pointed out in Ref. l27l. 
We expect that if the RG could be carried out exactly, 
we would find no divergent susceptibilities in this case. It 
would be very interesting to find systematic extensions of 
the approximate RG analysis used here which would be 
powerful enough to capture the Mcrmin- Wagner physics. 

We also considered the effects of an applied perpen- 
dicular magnetic field on the system when it is in the 
antiferromagnetic phase. In this case the fluctuations 
effects are weaker than at B = since the broken con- 
tinuous symmetry is the U(l) subgroup of the full spin 
SU(2) group. At B ^ a finite temperature transition 
into a power-law correlated state is in fact possible. Our 
variational mean field investigation was motivated by the 
fact that we find an antiferromagnetic state in our RG 
calculations for short-ranged interactions, as well as by 
experimental data^ on the gap size as a function of an 
applied magnetic field. We find that the antiferromag- 
netic order parameter increases quadratically with the 
field for low fields, then acquires a dependence of the form 
B/\n(B/Bo) for large fields. We also determined the 
gap by considering the energy required to create particle- 
hole excitations of our variational ground state. We find 
that this energy is the sum of the energies of the par- 
ticle and the hole given by the single-particle "auxiliary 
spectrum" , plus a term linear in the magnetic field. The 
excitation that gave us the smallest such energy was as- 



sumed to determine the energy gap in the system. We 
found that the gap has a slight non-monotonic behav- 
ior for low fields, followed by a quasi-linear increase at 
higher fields. We also compared this prediction to the 
experimental data and found that good agreement can 
be achieved. 

One reason for our switch to mean-field methods for 
treating this problem is that we have already established 
via RG methods the presence of the antiferromagnetic 
state under certain conditions. As long as we are consid- 
ering a case in which we know this phase to be present, 
and because said phase is gapped, we expect that an 
expansion around the mean-field solution will be conver- 
gent, thus justifying our use of such methods in studying 
the phenomenology of the phase. Another reason for our 
use of mean-field methods as opposed to the RG methods 
used previously in the zero-field case is the fact that the 
energy spectrum for the non-interacting problem is dis- 
crete, rather than continuous, and momentum k is not a 
good quantum number, making the use of RG methods 
more difficult. While we expect our mean- field methods 
to be fairly accurate, such methods are still only approx- 
imate. The problem of developing a more sophisticated 
technique for determining the AF order parameter and 
the energy gap in the system is therefore of great interest. 
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Appendix A: Long-distance behavior of the screened 
Coulomb-like interaction 

We will now demonstrate that Equation (f3T)| is the 
correct long-range behavior of the screened Coulomb-like 
interaction, Equation (|2"9"|) . We start by rewriting the 
sum using the identity, 

-i n />oo 

t = J-\ due~ r2u , (Al) 
r Jo 

obtaining 

V{B) = -==U Q l due- u2 W® 2 Y (-l) n e- n2u \ 

(A2) 

We may evaluate the sum in terms of the Jacobi theta 
function, 

oo 

d i {z,q)= i-l) n q n2 e 2mz , (A3) 



to obtain 

V(R) = -=U due-" 2(fl/ « )2 i9 4 (0,e- u2 ). (A4) 

V 71 " JO 

We now use the identity^!, 



p (4 Z 2 +7r 2 )/41og g y^ p fc(fc+l)7r 2 /log 9 

V- log g 

(2fc + 1)ttz 



fc=0 



x cosh 



logg 



(A5) 



so that 



V(R) =4C7 V / du-e- ul 



(A6) 

This integral can be evaluated in terms of modified Besscl 
functions of the second kind; the result is 



V{B) = AU Y j K 



k=0 



(2k + 1)tt 



(A7) 



For large values of x, the modified Bcssel function K n (x) 
can be approximated as 



K n (x) 



(A8) 



We see that, in the above sum, the most dominant term 
for R ^> £ is the k = term, since the values of suc- 
cessive terms decrease exponentially with increasing k. 
Therefore, we arrive at the form quoted in the main text, 
Equation ((30)) . 



Appendix B: Solution of the non-interacting 
Hamiltonian in a magnetic field with an AF order 
parameter 

We now provide the details of the diagonalization of 
Hq. To diagonalizc Ho, we first note that the source 
term is diagonal in layer, valley, and spin space, while 
the "kinetic" term is only diagonal in valley and spin 
space. This allows us to split the problem into the diago- 
nalization of four 2x2 matrices; our wave functions will 
have a definite valley pseudospin and real spin orienta- 
tion. Let us consider the block corresponing to the +K 
valley and spin up. We must solve 



A 

(7T x +i7r y ) 2 
2m* 



(7T X — ITTy) 

2m* 

-A 



i/>(x,y) =E1>(x,y), (Bl) 



where ip(x, y) is a two-component spinor corresponding 
to the +K valley and spin up components of the full 
eight-component spinor; the other six components are 
all zero. We will work in the Landau gauge, in which 
B = —BySc. For this gauge, the above becomes 



1 

2m* 



A 

(Px + ^-y 



eB r 



iPy 



-A 



IPy) 



ip(x,y) 



Eip(x,y). 

(B2) 



Let us now assume the following form for tp(x,y): 

1 



ip(x,y) 



Akx 



®k(y), 



(B3) 



where $fe(y) is another two-component spinor. Upon 
substitution into our equation, we obtain 



A 



^(k+fy + ip y y 



^(k+zfy-iPy) 



E$ k (y). 

(B4) 



We may now write this in terms of the operators, 

k 



+ y + %- 



Py 
m*ui c 



(B5) 



-^f- is the cyclotron frequency of the elec- 



where w c 

trons in our system. These operators may be verified to 
satisfy the commutation relation, [afc,a£] = 1. In terms 
of these operators, the equation becomes 

A 



u c (a\) 



^2 

A 



$ fc (y) = E$ k (y). 



(B6) 



Let us now define normalized functions 4>k,n(y) such 

that a\ak4>k,n(y) = n(f>k,n(y)> these will just be the 
usual harmonic oscillator-like wave functions that emerge 
in the solution of the free electron gas in a magnetic 
field. These functions may be shown to satisfy the 
r elation s, a k <j) k>n (y) = y/ncf) ktn -i(y) and a^ fei „(y) = 
yn + 10fc n+i(y)- We now assume the following form for 



$k(y) 



ak,n4>k,n(y) 
Pk,n4>k,n-2(y) 



(B7) 



We find that this form satisfies our equation, provided 
that a k ,n and j3 k _ n satisfy 



A (xi c yjn(n — 1) 

ijj c \/ n(n — 1) —A 



ctk, 



= E 



Otk,n 
(3k,n 



(B8) 

and |afc, n | 2 + |/3fc, n | 2 = 1- We have thus reduced the 
problem to solving for the eigenvalues and eigenvectors 
of a 2 x 2 matrix. The eigenvalues arc E = ±E n , where 



E n = v/n(n-lK 2 + A 2 , (B9) 
and the corresponding eigenvectors are given by 
1 



(BIO) 



In the — K valley, the positions of the creation and anni- 
hilation operators a\ and au will be interchanged; in this 
case, we must instead assume that 



<Xk,n<l>k,n-2(y) 
(3k,n4>k,n(y) 



(Bll) 



This will give us the same eigenvalue problem as before. 
For the spin down case, we simply reverse the sign on A in 
our equations; we obtain the same eigenvalues as before, 
but and will switch values. This implies that, 
at least for n > 2, each of our energy levels is four- fold 
degenerate, due to valley and spin degeneracies. For n = 
or 1, on the other hand, there is no valley degeneracy 
for a given spin. In these cases, the eigenfunctions are 



$fc(y) 



in the +K valley and 



$k(y) 











(B12) 



(B13) 



in the — K valley. The former corresponds to the energy 
eigenvalue, E = A, while the other corresponds to E = 
—A. Each of these levels is still four-fold degenerate, but 
this time due to orbital and spin degeneracies. Note that 
we may still use the wave functions quoted earlier even 
for this case if we adopt the convention that 4>k.n{y) is 
identically zero if n < 0. 



Appendix C: Derivation of the variational mean-field 
equation for an AF order parameter 

We now describe how to derive the variational mean- 
field equation, Eq. (|4Tj) . We first provide details on the 
calculation and minimization of the ground-state expec- 
tation value of the interacting Hamiltonian with an anti- 
ferromagnetic order parameter. Throughout this deriva- 
tion, we will assume that A > 0. We start by rewriting 
the field operators in our Hamiltonian in terms of the 
eigenstates of the Hq derived above; the formula for this 
is given by Equation (f3"T| . If we label the positive energy 
eigenvalues as 2?k, niT)S = E n = \Jn{n — l)w 2 + A 2 , then 
Hq becomes 

#0 = ^ Ek,n,T,s( a k,n,T,s a k,n,T,s + b l,n,T,s b k,n,T,s) 
k,n,r,s 

— 2^1 Ek,n,T,s- (CI) 

The expectation value of Hq with respect to the ground 
state |0) is then 

(0|Fo|0> = - E k , n , T , s = -4dJ2E n -4d\A\, 

k,n,T,s n 

(C2) 



where d is the degeneracy of each Landau level due to 
the wave number k. 

Now we turn our attention to the interaction, starting 
with the quadratic term, 



A / d 2 r^(r)l 2 a z s z tp(r). (C3) 



H { P 

In terms of the eigenstates of Hq, this becomes 

H? ] = -A j rf 2 r^i(r)l 2ff 2 S z i(r) fl La n . (C4) 

ran 

Here, we let the indices m and n stand for all of the 
quantum numbers characterizing a given eigenstate, and 
we use a m to stand for either a positive or negative energy 
state, with the understanding that the negative energy 
state is given by a m = b^. Upon taking the expectation 
value of this term with respect to the ground state, we 
find that only the negative energy states contribute: 

(0\HP |0) = -A f d 2 v J2^ m )Hr)h<y z s z yj m (r) 

J m 

(C5) 

Evaluating all of the sums and integrals using the expres- 
sions given above for the wave functions, we may write 
this as 

(0| H [ p |0) = -i,L a ATr(l 2 a z S z S_), (C6) 

where 

j n 

= Tlrr^ciNls + t z <t z 1 2 - t z 1 2 <j z - Y1 2 cj z s{,§7) 
4ir£ B 

Evaluating the trace, we obtain 



(C8) 



where Ib = \Jc/eB is the magnetic length, L x and L y 
arc the dimensions of our system, and 



n>2 



A_ 



1. 



(C9) 



We now consider the quartic terms. Each of these 
terms, neglecting the coupling constants and integrals 
over position, has the form, 



(CIO) 



[^(r)S^(rr, 



where S is a matrix. Substituting in Equation (|37p , and 
adopting the same conventions as before, this becomes 

£ [^l(r)SMr)][^(r)S^(r)]ala n ala q . (Cll) 



We now take the expectation value of this expression with 
respect to the ground state. This expectation value will 



involve the expression, (0| aj„a„a^a 9 |0). The only way 
for this to be non-zero is if m and q are negative-energy 
states. We must also require that n and p either be both 
positive-energy states or both negative-energy states. In 
the former case, we obtain, using the anticommutation 
relations for fermions, 

(0| b m a n a p b\ |0) = 8 mq 6 np . 

In the latter case, we obtain 

(0| b m b\b p b\ |0) = S mn S pq . 

Putting these results together, the above quartic form 
becomes 

E l(i> m )Hr)Sri (r)] [(V>+ y (r)Sty- (r)] 

m,n 

i 



+ E[(V'™) t (r)5 , V'™(r)][(V p -) t (r)^-(r)], (C12) 

Upon evaluating all sums and integrals, this becomes 

Tr(S*S + S'S_) + [Tr(S'S_)] 2 , (C13) 

where 

m 

= T-^-^o(Nl 8 +T z a z l2 + r z l2a z + Yl 2 a z ^.4) 

Evaluating the traces, we find that the total energy 
E V!l r = (Q\H\0) of our system is, noting that d = 
L x L y /(2Ttl%), 



2L x L y 



E E » + A 



\L x L y g Alg [ —jq- ) , 



2N 



2L L 

— |pAF + 4L x L y (g Alg + g A2n + Ag El< ] 



N 
4< 



/ Y 



(C15) 



where 



n>2 



(C16) 



We will now minimize this energy with respect to A. First, we take the derivative of the above expression, which may 
be written as 



dE v 



dA 



2L X Ly 



A 



1 



Anil 



i.9A lg +9A 2u +AgE K )Y 



jr.. 



n>2 



1 - 



E 2 



(C17) 



We now set this derivative to zero. We note that the 
second factor can never be zero, since E n > A for all n > 
2, and therefore it is always positive. We may therefore 
drop this factor. The remainder, upon simplifying and 
noting that Ib = 1 / y/m*u) c , thus yields Eq. (|38j) . 

We now wish to rewrite Eq. (|38[) in such a way that 
we may send the upper cutoff to infinity. We start be 
rewriting the sum as an integral over a "Dirac comb" : 



A 



u c / dX V S(X - n) 

= (C18) 

ffeff 



We now use the identity, 

oo oo 

E *(A-n)= E 



(C19) 



to obtain 



dA ■ 



1 



+2w c E / 

, J3 



3/2 ^X(X - T)U* + A 2 

' cos(27rfcA) 



^J 3 /2 v/A(A - 1)uj 2 + A 



■ A" 
1 



5off 



(C20) 



fc— — oo 



The integral that appears in the sum converges for all 
ft > 1, so we may already send the upper limit to in- 
finity. However, we must still take care of the first inte- 
gral, which will diverge if we do the same with it. We 
may use the zero- field equation, Equation (|40[) , to rewrite 
the right-hand side such that we eliminate g e s from the 
equation. If we do this and also introduce the change 
of variables, e 2 = A(A — 1)lu 2 , into the first integral, we 
may then move the first integral to the right-hand side, 
obtaining an integral that converges if wc send the upper 



limit to infinity, namely 
1 



uj c I de 
lo 



+UJ C 



V^Ta 2 ^ e 2 + 1^2 V^T 

V3/2 e 1 



de 



e 2 + l w a Ve 2 + A2 



(C21) 



cut. The contribution from the circular arc will vanish 
as we increase the radius to infinity since the integrand 
decreases exponentially as we do so. This leaves only 
contributions from the radii. The contribution from the 
radius along the real axis is just the integral that appears 
in the equation. This means that the contribution from 
the radius parallel to the imaginary axis is equal to this 
integral. We may therefore write 



where Ao is the value of A at zero magnetic field. Upon 
evaluating this integral, our equation becomes 

oo^ C0S (27r/cA) u c 

2uj c > / dX y ' + -2 

y/\(\-l)u* + A* A 



/OO 
(l.T 



: 2 + a 2 - ± 



Rc 



= In 



A 



A 2 + \u>l 
A 



g — 2izkx r 

i I dx' 

° J(l + ix') 2 + a 2 - \ 



This equation may be rewritten in terms of the dimen- 
sionless parameters, a = A/uj c and /? = cj c /Aq: 



(C22) Note that we dropped a factor of e~ 27 " fe ; since k is an 
integer, this factor is always equal to 1. If we substitute 
this back into the equation, we find that the sum on k is 
just a geometric series with a common ratio of ~ e ~ 2nx . 
We may therefore perform the summation, obtaining 



lV / dX 

= ln/3. 



cos(2irkX) 



I -In 
a 



a z + 



1(a) + 



1 



In 1 



ln/3, 



(C24) 



(C23) 



We can see that the left-hand side is a monotonically- 
decreasing function of a for a > 0. In fact, as a — > + , 
the left-hand side increases indefinitely due to the second 
term, while, as a — > oo, the expression decreases indefi- 
nitely due to the third term. This equation therefore has 
a single positive solution for a for any given value of p. 
While we would need to solve the equation numerically 
for general values of j3, we can derive approximate solu- 
tions analytically for very large and very small values of 
/?. 

Before we do this, however, let us first rewrite this 
equation in an equivalent form. We start by changing 
variables in the first term to x = X — i , obtaining 



where 



1(a) = 2 Re 



dx' 



o ^(1 + ix') 2 + a 2 



l e 

4 



27TX' 



(C25) 



2^(-l) fc Rc / dx 
i — i J l 



We can see that the integral 1(a) converges for all values 
of a; the integrand is analytic everywhere on the interval 
of integration and decreases exponentially for large x' . 
At large values of a, we can show that this integral falls 
off as a~ 3 . We first note that the integral is dominated 
by small values of x' due to the Fermi occupation factor- 
like expression. With this in mind, we may pull out a 
factor of a from the square root, obtaining 



J2nikx 



k=l 

In | 1 



1 

a 



-2Rc 



dx' 



1 



(1 + ix') 



In -1/2 
4 



1 



i 



ln/3. 



We now note that, as a function of x, the integral is 
analytic in the entire complex plane, except for a branch 
cut. For a < i this branch cut can be chosen to be 



on the real axis and in the interval 



a 2 < x < 



Since a is large, we now have a small parameter with 
respect to which we may perform an expansion of the 
square root. The constant term in this expansion gives no 
contribution, since the total result will be purely imag- 
inary. The lowest-order non-zero contribution will, in 
fact, be given by 



-j — a 2 . For a > ^, on the other hand, the branch 

cut may be chosen to lie along the imaginary axis. In 
either case, we may integrate this function over a large 
quarter circle centered at the point, x = 1, in the complex 
plane and with one of the radii along the positive real 
axis and the other parallel to the positive imaginary axis 
and obtain zero since we will always avoid the branch 



dx' 



1 



-,2-nx' 



1 



24a 3 ' 



We see that this term is of the order a -3 , as asserted 
earlier. 

We may derive a good closed-form approximation to 
1(a) as follows. Let us first expand the square root in the 



integrand in powers of x' . To the lowest no n- vanishing 
order, we obtain 



/>oo 

1(a) « -2 / dx' 
Jo 



1 



1 + a 2) 3/2 e 2 ™' + I 3(3 + 4a 2 )3/2 



We now rewrite this expression so that its value at a = 
matches the exact value of 1(0). Doing so, we obtain 



1(a) « - 



[( l/(0))-2/3 + 4 . 62/3 a 2]3/2 ■ 



(C26) 



If we were to plot this expression alongside the exact 
expression for 1(a), then we would see that it follows 
the exact expression very closely. In fact, if we use this 
expression to solve Equation (f4"Tj) . then the solution that 
we obtain is very close to the solution obtained from the 
exact 1(a). 



Appendix D: Derivation of the excitation spectrum 

We now present the details of our derivation of the ex- 
citation spectrum, given by Eq. (j4"T|). As stated before, 
we begin by constructing a particle-hole excitation of our 
trial ground state, b^a^ |0). We then find the difference in 
the expectation value of the Hamiltonian, given by Eq. 
(|32|) . between the excited state and the ground state; 
this is taken to be the excitation energy Throughout 



this calculation, we assume that the AF order parame- 
ter A > 0. Let us begin with the quartic terms. If we 
take the difference in expectation value of these terms 
between the excited state and the ground state, we find, 
after straightforward but tedious application of anticom- 
mutation relations and dropping terms that will vanish 
in the thermodynamic limit, that the contribution to the 
excitation energy is 

5E A = J2 <?s[Tr(S£ Q/3 )Tr(SI]_) + ±Tr(S£ + _S£ Q(3 )], 
s 

(Dl) 

where 

£+- = E^n( r )^n) t ( r )-^( r )^n) t W]' ^ 
n 

V a p = J d 2 v[^(v)(^)Hv)~i,- p (v)(i,- p )\v)]. 

(D3) 

We may now evaluate the sums and integrals in the above 
expressions, obtaining 



£+- = ^-^c(Yl 2 <J z s z + t z 1 2 s z ) (D4) 

Z7T 



and 



( 1 - si^— ) (1 + nr a )(l - (7*)(1 + sis z) o 



( A \ , ( A 

= TE \l + S 1 —\(l + TlT z )(l + a z )(l + SlS z )e ni -l +Tl , 

- i f 1 - S2 A) (i + T2 t 2 )(1 + <7 Z )(1 + S2S z )6 nz -i +T2 - i (l + s 2 — ) (1 + t 2 t z )(1 - a z )(l + s 2 s z )0 n2 -i- T2 



(D5) 



Here, 6 n is if n < 0, or 1 otherwise. We may now 
evaluate the traces in Eq. (|D1[) . Doing so, and using the 
fact that 

7Tb 

^c(gA lg + g A - 2u + 4g EK )Y = A, (D6) 

we find that 

( E ni E n2 ) 
+ -^ u) c(9A lg + gA 2u - 4 g£;jc )(TiSi - r 2 s 2 ). 

(D7) 

We now consider the quadratic terms. One of these terms 
simply gives us the single-particle "auxiliary spectrum" , 



while the other is the quadratic term that we chose to 
consider as part of the interaction term. We find, upon 
application of anticommutation relations as before, that 
the contribution from these terms to the excitation en- 
ergy is 

5E 2 = E ni + E n2 - ATr(l ( T 3 s z Sa/j). (D8) 
Upon evaluating the trace, we obtain 

5E 2 = E ni + E n2 - A 2 f + . (D9) 

Combining these two contributions, we arrive at Eq. (|47|) . 
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